In [1]:
from itertools import izip
from time import time
import numpy as np
import astropy
from pearce.mocks.customHODModels import *
from pearce.mocks import cat_dict
from scipy.optimize import minimize
In [2]:
from SloppyJoes import lazy_wrapper
In [3]:
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
In [4]:
AB = True
In [5]:
PRIORS = {'f_c': (0, 1),
'alpha': (0, 2),
'logMmin':(10,14),
'logM1': (10, 15),
'logM0': (9,15),
'sigma_logM': (0.3, 1.5),
'logMcut': (9,15),
'logMlin':(9,15),
'f_cen': (0.0,1.0)}
_cens_model = RedMagicCens
cens_model = _cens_model(z = 0.0)
#cens_model = AssembiasReddick14Cens()
_sats_model = RedMagicSats
#sats_model = AssembiasReddick14Sats()
cosmo_params = {'simname':'chinchilla', 'Lbox':400.0, 'scale_factors':[0.658, 1.0]}
cat = cat_dict[cosmo_params['simname']](**cosmo_params)#construct the specified catalog!
cat.load(1.0, HOD=(_cens_model, _sats_model), hod_kwargs = {'cenocc_model': cens_model})
LBOX = 400.0
#sats_model.modulate_with_cenocc = False
In [6]:
cat.model.model_dictionary
Out[6]:
In [7]:
cens_model = cat.model.model_dictionary['centrals_occupation']
sats_model = cat.model.model_dictionary['satellites_occupation']
In [8]:
def resids(theta,params,cens_occ, sats_occ,mbc):
cens_model.param_dict['f_c'] = 1.0
sats_model.param_dict['f_c'] = 1.0
cat.model.param_dict['f_c'] = 1.0
cens_model.param_dict.update({p:x for p, x in izip(params, theta)})
sats_model.param_dict.update({p:x for p, x in izip(params, theta)})
cat.model.param_dict.update({p:x for p, x in izip(params, theta)})
cens_preds = cens_model.mean_occupation(prim_haloprop = mbc)
sats_preds = sats_model.mean_occupation(prim_haloprop = mbc)
#Weird edge cases can occur?
cens_preds[cens_preds < 1e-9] = 0
sats_preds[sats_preds < 1e-9] = 0
cens_vars = cens_preds*(1-cens_preds)+1e-6
sats_vars = sats_preds + 1e-6
Ngal_pred = np.sum(cens_preds+sats_preds)
Ngal_obs = np.sum(cens_occ+sats_occ)
idx = sats_occ > 0
#log_sats_diff = (np.log10(sats_preds) - np.log10(sats_occ) )
#log_sats_diff[np.isnan(log_sats_diff)] = 0.0
#log_sats_diff[log_sats_diff == -np.inf] = 0.0
#log_sats_diff[log_sats_diff == np.inf] = 0.0
return np.r_[ (cens_preds-cens_occ),sats_preds-sats_occ, np.array([Ngal_pred-Ngal_obs]) ]
#return np.r_[cens_preds[0,:]-cens_occs[0,:], Ngal_pred-Ngal_obs]
In [9]:
catalog = astropy.table.Table.read('/u/ki/swmclau2/des/AB_tests/abmatched_halos.hdf5', format = 'hdf5')
In [10]:
mag_cut = -21
min_ptcl = 200
if AB:
catalog = catalog[np.logical_and(catalog['halo_mvir'] > min_ptcl*cat.pmass, catalog['halo_vpeak_mag'] <=mag_cut)]
else:
catalog = catalog[np.logical_and(catalog['halo_mvir'] > min_ptcl*cat.pmass, catalog['halo_vvir_mag'] <=mag_cut)]
In [11]:
if not AB:
pass
else:
MAP = np.array([ 12.87956269, 12.24461447, 0.5345765, 13.98105124, 1.04527197])
['$\\log{M_{min}}$', '$\\log{M_0}$', '$\\sigma_{log{M}}$', '$\\log{M_1}$', '$\\alpha$']
names = ['logMmin', 'logM0', 'sigma_logM', 'logM1', 'alpha']
hod_params = dict(zip(names, MAP))
In [12]:
ab_params = {'mean_occupation_centrals_assembias_param1':0.4, 'mean_occupation_satellites_assembias_slope1':3,\
'mean_occupation_satellites_assembias_param1':-0.5, 'mean_occupation_centrals_assembias_slope1':3,}
In [13]:
sats_model.param_dict.update(cens_model.param_dict)
In [14]:
param_dict = hod_params
#param_dict.update(ab_params)
cens_model.param_dict.update(param_dict)
sats_model.param_dict.update(param_dict)
params = sats_model.param_dict.keys()
########################
params.remove('f_c')
#######################3
ndim = len(params)
In [15]:
halo_table = cat.halocat.halo_table[cat.halocat.halo_table['halo_mvir'] > min_ptcl*cat.pmass]
In [16]:
detected_central_ids = set(catalog[catalog['halo_upid']==-1]['halo_id'])
In [17]:
from collections import Counter
def compute_occupations(halo_table):
#halo_table = cat.halocat.halo_table[cat.halocat.halo_table['halo_mvir'] > min_ptcl*cat.pmass]
cens_occ = np.zeros((np.sum(halo_table['halo_upid'] == -1),))
#cens_occ = np.zeros((len(halo_table),))
sats_occ = np.zeros_like(cens_occ)
detected_central_ids = set(catalog[catalog['halo_upid']==-1]['halo_id'])
detected_satellite_upids = Counter(catalog[catalog['halo_upid']!=-1]['halo_upid'])
for idx, row in enumerate(halo_table[halo_table['halo_upid'] == -1]):
cens_occ[idx] = 1.0 if row['halo_id'] in detected_central_ids else 0.0
sats_occ[idx]+= detected_satellite_upids[row['halo_id']]
return cens_occ, sats_occ
In [18]:
from halotools.utils.table_utils import compute_prim_haloprop_bins
def compute_hod(masses, centrals, satellites, mass_bins):
mass_bin_idxs = compute_prim_haloprop_bins(prim_haloprop_bin_boundaries=mass_bins, prim_haloprop = masses)
mass_bin_nos = set(mass_bin_idxs)
cens_occ = np.zeros((mass_bins.shape[0]-1,))
sats_occ = np.zeros_like(cens_occ)
for mb in mass_bin_nos:
indices_of_mb = np.where(mass_bin_idxs == mb)[0]
denom = len(indices_of_mb)
#TODO what to do about bout 0 mean std's?
cens_occ[mb-1] = np.mean(centrals[indices_of_mb])
sats_occ[mb-1] = np.mean(satellites[indices_of_mb])
return cens_occ, sats_occ
In [19]:
mass_bin_range = (9,16)
mass_bin_size = 0.1
mass_bins = np.logspace(mass_bin_range[0], mass_bin_range[1], int( (mass_bin_range[1]-mass_bin_range[0])/mass_bin_size )+1 )
mbc = (mass_bins[1:]+mass_bins[:-1])/2
In [20]:
cens_occ, sats_occ = compute_occupations(halo_table )
mock_masses = halo_table[halo_table['halo_upid']==-1]['halo_mvir']
#mock_concentrations = halo_table[halo_table['halo_upid']==-1]['halo_nfw_conc']
In [21]:
cen_hod, sat_hod = compute_hod(mock_masses, cens_occ, sats_occ, mass_bins)
In [22]:
param_dict.keys()
Out[22]:
In [23]:
params
Out[23]:
In [24]:
vals = np.array([param_dict[key] for key in params])
cens_idxs = halo_table['halo_upid'] == -1
args = (params, cen_hod, sat_hod,mbc)
print params
In [25]:
test = cens_model.mean_occupation(prim_haloprop = cat.halocat.halo_table['halo_mvir'][:100],\
sec_haloprop= cat.halocat.halo_table['halo_nfw_conc'][:100])
print np.mean(test)
In [26]:
mbc.shape
Out[26]:
In [27]:
resids(vals, *args)
Out[27]:
In [28]:
lazy_wrapper(resids, vals, func_args = args,maxfev = 500, print_level = 1, artol = 1e-6)
Out[28]:
In [29]:
print params
In [30]:
print MAP
In [ ]: